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This  paper  examines  the  effect  of  the  electrical  double  layer  on  the  performance  of  a  lithium  ion  battery 
electrochemical  cell.  We  begin  by  introducing  the  Poisson  Nernst— Planck  equations  of  electrochemistry 
to  describe  ion  transport  within  a  representative  liquid  solvent  and  derive  an  expression  for  the  current 
—voltage  relationship  in  the  electroneutral  liquid  within  the  separator  pores.  Different  assumptions 
about  the  electrical  double  layer  lead  to  variation  of  the  lithium  ion  concentration  profiles  in  the  liquid 
electrolyte,  which  alter  the  cell  voltage  during  discharge.  The  contribution  of  the  electrical  double  layer 
to  the  cell  overpotential  is  combined  with  the  bulk  liquid  potential  difference  and  a  simplified  treatment 
of  the  electrode  solid  phase  to  obtain  an  expression  for  the  time-varying  cell  terminal  voltage.  We 
conclude  by  presenting  experimental  data  for  a  cell  using  a  graphite  anode  and  lithium  iron  phosphate 
cathode  to  validate  the  new  model,  which  includes  electrical  double  layer  effects,  and  construct  modi¬ 
fications  to  the  model  to  compensate  for  resistive  effects  that  are  associated  with  the  cathode  solid 
phase. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Modeling  is  extremely  important  to  advancing  the  science  of 
electrochemical  energy  storage.  Direct  and  dynamic  measurements 
of  concentration  and  electric  potential  profiles  within  a  Lithium-ion 
(Li-ion)  cell  are  currently  not  possible  due  to  the  micro-scale 
physical  dimensions  of  the  battery  electrodes  and  separator.  As  a 
result,  researchers  have  relied  on  modeling  as  a  means  of  under¬ 
standing  the  complex  processes  governing  electrochemical  cells. 
There  have  been  a  number  of  fundamental  studies  of  Li-ion  battery 
operation  based  on  porous  electrode  theory  [1]  since  the  original 
work  modeling  a  lithium-polymer  insertion  cell  [2],  Porous  elec¬ 
trode  theory  is  an  extremely  powerful  methodology  which  ac¬ 
counts  for  the  presence  of  distinct  phases  within  the  battery  using 
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1  Units  within  the  Nomenclature  section  denote  the  units  used  for  dimensional 
quantities  indicated  by  the  (*)  superscript.  If  the  (*)  superscript  is  absent  then  the 
variables  are  dimensionless. 
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superposition.  Detailed  knowledge  of  the  surface  morphology  of 
electrodes  and  separator,  which  is  difficult  to  obtain  and  compu¬ 
tationally  prohibitive  to  model,  is  not  required.  The  solid  micro¬ 
structure  of  both  electrodes  and  separator  is  accounted  for  only 
with  a  pore  volume  fraction  and  Bruggeman  coefficient  for  tortu¬ 
osity.  Therefore,  concentration  and  potential  profiles  are  computed 
in  a  volume-averaged  sense,  and  models  based  on  porous  electrode 
theory  present  a  macroscopic  interpretation  of  battery  electro¬ 
chemistry.  The  theory  is  flexible  enough  to  allow  for  the  inclusion 
of  features  that  are  unique  to  a  given  chemistry,  such  as  phase 
change  during  intercalation,  path  dependence,  or  resistive-reactant 
nature  of  electrode  active  materials.  All  of  these  phenomena  have 
been  postulated  to  occur  in  the  graphite/lithium  iron  phosphate 
couple  that  is  studied  in  this  paper  [3-6]. 

Since  porous  electrode  models  compute  averaged  quantities 
over  a  region  that  is  small  with  respect  to  the  overall  electrode 
dimensions  but  large  compared  to  the  pore  structure  [1],  Li-ion 
battery  modeling  literature  typically  applies  the  assumption  of 
electroneutrality  within  commonly  used  non-aqueous  liquid  elec¬ 
trolytes.  This  assumption  is  justified  based  on  the  assertion  that 
charge  separation  over  a  macroscopically  significant  distance  (the 
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Nomenclature1 

14 

liquid  phase  voltage  (V) 

Vrr 

resistive  reactant  overpotential  (V) 

a 

Butler-Volmer  transfer  coefficient 

X 

separation  of  variables  function  of  space 

y 

active  material  volume  fraction 

X 

spatial  variable  in  direction  across  unit  cell  thickness 

e 

ratio  of  Debye  length  to  separator  thickness,  X/Lsep 

(m) 

ee 

electrical  permittivity  (F  m  ’) 

y 

inner  solution  spatial  variable,  x/e 

? 

zeta  potential,  or  equivalently  the  potential  difference 
across  the  diffuse  layer  (V) 

z 

ion  valence 

X 

Debye  length  (m) 

Subscripts 

electric  potential  (V) 

00 

denotes  a  property  at  the  interface  between 

A 

current  collector  area  (m2) 

electroneutral  liquid  and  EDL 

C 

current  magnitude  in  terms  of  cell  capacity  (1 C  =  2.3  A) 

A 

refers  to  a  cation  (Li+)  property 

C 

concentration  (mol  m  3) 

B 

refers  to  an  anion  (PFg )  property 

Co 

electroneutral  concentration  (mol  m“3) 

a 

C/t,an/C/\,an,max 

D 

diffusion  coefficient  (m2  s-1) 

an 

denotes  a  property  of  the  anode 

F 

Faraday  constant  (C  mor1) 

b 

CA.cJCA,a.max 

I 

external  current  demand  (A) 

ca 

denotes  a  property  of  the  cathode 

J 

current  density  (A  itt2) 

EDL 

property  of  the  electrical  double  layer  in  its  entirety, 

Jd 

diffusion  limited  current  density  (A  m-2) 

including  Stem  and  diffuse  layers 

Js 

solid  boundary  flux  scale  (A  m-2) 

i 

refers  to  a  cation  property  if  i  =  A  or  an  anion  property 

chemical  reaction  rate  constant  for  oxidation  reaction 

if  f  =  B 

(A  m  mor1  s_1) 

j 

refers  to  an  anode  property  if  j  =  an  or  a  cathode 

kr 

chemical  reaction  rate  constant  for  reduction  reaction 

property  if  j  =  ca 

(A  m4  mol-2  s_1) 

max 

denotes  the  saturation  concentration  value 

L 

thickness  of  separator,  anode,  or  cathode  (m) 

s 

Stern  layer  property 

n 

Q 

index  for  inifinite  series  solution 
capacity  removed  (Ah) 

sep 

refers  to  a  separator  property 

R 

universal  gas  constant  (J  mol-1  K) 

Superscripts 

Rc 

contact  resistance  (Q  m2) 

+ 

denotes  an  ion  with  valence  equal  to  1 

Rrr 

resistive  reactant  resistance  (Q  m2) 

- 

denotes  an  ion  with  valence  equal  to  -1 

SOC 

state  of  charge 

dimensional  quantity 

T 

temperature  (K) 

~ 

refers  to  the  non-equilibrium  portion  of  EDL  electric 

T 

separation  of  variables  function  of  time 

potential 

t 

time  (s) 

refers  to  electric  potential  defined  relative  to  the 

u 

V 

open-circuit  voltage  relative  to  Li/Li+  (V) 
terminal  voltage  (V) 

adjacent  bulk  potential 

volume  averaging  distance)  would  require  a  prohibitively  large 
electric  field  [7],  However  at  the  boundary  between  the  liquid  and 
solid  phases,  an  interfacial  region  known  as  the  electrical  double 
layer  (EDL)  exists  where  the  assumption  of  electroneutrality  no 
longer  holds  [8],  A  method  of  accounting  for  capacitive  effects  of 
the  EDL  in  cell  terminal  voltage  has  previously  been  presented  [9], 
The  scope  of  Ref.  [9]  did  not  encompass  some  important  features  of 
the  EDL  such  as  the  coupling  between  ion  concentration  and  po¬ 
tential  described  in  several  classical  works  [10—13]  Given  the 
general  lack  of  literature  pertaining  to  the  EDL  as  it  applies  to  Li-ion 
batteries,  we  feel  it  is  appropriate  to  examine  this  feature  in  greater 
detail. 

A  more  accurate  understanding  of  the  EDL  is  important  for  its 
impact  on  interfacial  charge  transfer.  Porous  electrode  models 
postulate  that  charge  transfer  kinetics  are  governed  by  the  poten¬ 
tial  difference  between  the  solid  phase  and  the  electroneutral 
portion  of  the  liquid  outside  the  EDL  [7],  This  assumption  is  made 
implicitly  in  the  application  of  Butler— Volmer  kinetics,  which  is 
treated  as  a  semi-empirical  relationship  (though  possessing 
fundamental  origins)  between  the  local  intercalation  current  and 
the  potential  difference  between  solid  and  liquid  relative  to  the 
open-circuit  voltage  versus  a  hypothetical  lithium  reference  elec¬ 
trode.  However,  the  Stern  layer  is  the  location  of  the  reaction  plane 
and  the  potential  difference  across  it  represents  the  actual  activa¬ 
tion  energy  barrier  for  the  intercalation  reaction  [10],  Thus  it  is  of 


fundamental  interest  to  consider  only  the  Stern  layer  potential 
difference  when  examining  the  current-voltage  relationship  of 
electrode  kinetics.  The  two  cases  of  EDL  physics  are  that  the  entire 
potential  difference  occurs  within  either  the  Stern  layer  or  the 
diffuse  layer.  A  complete  methodology  for  dividing  the  potential 
difference  of  the  EDL  between  the  diffuse  layer  and  the  Stern  layer 
is  discussed  in  Ref.  [10]  in  relation  to  a  micro-battery,  and  for 
galvanic  cells  in  Ref.  [  14],  Results  are  presented  for  a  range  of  kinetic 
rate  constants  selected  to  examine  various  limits  of  battery  oper¬ 
ation  with  no  particular  system  studied.  For  an  intercalation  bat¬ 
tery,  the  effective  rate  constant  changes  over  the  course  of 
discharge  due  to  the  changes  in  active  material  lithium  content. 
This  phenomena  is  governed  by  the  composition  ranges  of  the 
electrodes  chosen  by  the  manufacturer  and  the  cell  state  of  charge 
(SOC),  and  provides  a  physical  mechanism  for  alteration  of  the 
effective  rate  constants  over  the  course  of  discharge. 

The  objective  in  this  paper  is  to  develop  an  analytical  model  that 
can  predict  battery  voltage  during  discharge.  The  presented  model 
structure  accounts  for  the  potential  difference  in  the  electroneutral 
liquid  within  the  separator  pores,  as  well  as  a  microscopic  inter¬ 
pretation  of  the  EDL,  which  leads  to  a  novel  view  of  electrode  ki¬ 
netics.  We  first  present  the  model  development,  where  we  make 
use  of  the  Poisson  Nernst— Planck  equations  to  model  liquid  phase 
transport  and  solve  a  representative  solid  diffusion  problem  in  each 
electrode.  Two  alternate  descriptions  of  the  EDL  potential  are  used 
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to  examine  the  impact  of  the  EDL  potential  on  the  cell  over¬ 
potential.  After  obtaining  analytical  solutions  for  all  variables  of 
interest  we  discuss  the  experimental  procedure  for  gathering  data 
from  commercially  available  cells  for  model  validation.  The  case 
study  uses  a  graphite  anode  and  iron  phosphate  cathode,  and 
electrolyte  composed  of  1  M  Li  PFg  dissolved  in  a  1:1  mixture  of 
ethylene  carbonate  and  dimethyl  carbonate.  We  compare  the 
model  results  with  galvanostatic  discharge  data  and  discuss  the 
areas  impacted  by  the  microscopic  hypotheses  of  the  EDL.  We 
introduce  additional  dynamics  that  are  proposed  to  be  related  to 
the  resistive  phenomena  within  the  cathode,  and  finally  present 
conclusions  and  areas  for  future  work. 


2.  Model  development 

To  develop  a  model  suitable  for  incorporating  the  effects  of  the 
EDL  on  the  battery  voltage  during  galvanostatic  operation,  we 
consider  the  transport  of  ions  that  occurs  along  a  single  liquid-filled 
pore  of  the  separator  that  terminates  at  solid  faces  of  the  anode  and 
cathode.  Fig.  1  shows  the  modeled  domain  of  this  work  in  the 
context  of  the  full  unit  cell  typically  modeled  in  literature. 


2.J.  Separator  liquid  phase 


The  assumptions  used  to  derive  the  governing  equations  are  as 
follows.  First,  significant  concentration  and  potential  gradients 
occur  only  in  the  direction  through  the  thickness  of  a  unit  cell, 
indicated  by  x  in  Fig.  1.  Second,  bulk  convection  of  the  electrolyte 
solvent  is  neglected.  Third,  the  external  current  is  independent  of 
time.  Finally,  the  liquid  phase  ionic  transport  occurs  via  steady  state 
conditions  since  the  time  scale  of  liquid  transport  is  much  less  than 
the  duration  of  the  capacity  test. 

The  Poisson  Nernst— Planck  (PNP)  equations  are  used  to  model 
transport  of  charged  species  in  dilute  solutions.  The  simplicity  that 
dilute  solution  theory  offers  compared  with  concentrated  solution 
theory  is  a  benefit  when  one  is  more  interested  in  factors  such  as 
the  descriptions  of  the  EDL  presented  in  the  introduction.  In  non- 
dimensional  form,  the  simplified  equations  that  result  from  the 
previously  listed  assumptions  are 


*ci  (dcid* 

9x2  dx  +  L'QX2) 


(1) 


where  c,-  gives  the  ionic  concentration  with  i  =  A  corresponding  to 
Li+  ions  and  i  =  B  denoting  PFg  ions,  z  is  the  valence  of  the  ion 
indicated  via  subscripts,  and  ip  is  the  electric  potential  in  the  liquid 
phase.  Conceptually,  Eq.  (1)  states  that  the  sum  of  fluxes  due  to 
diffusion  and  migration  are  balanced  at  steady  state.  The  electric 
potential  is  governed  by  Poisson’s  equation 


Cathode 
LibFeP04 
fVctive  Material 

°§ 

°o 

°o 

°o 

oR 


Active  Material 


e1^  =  ~(ZACA+ZBCB)  (2) 

The  parameter  e  =  A/Lsep  arises  from  converting  Eq.  (2)  to 
dimensionless  form.  The  parameter  A  is  defined  as 


and  called  the  Debye  length.  The  parameter  ee  is  the  permittivity  of 
the  liquid  solvent  taken  to  be  an  average  of  the  two  main  solvent 
components,  ethylene  carbonate  (EC)  and  dimethyl  carbonate 
(DMC).  The  Debye  length  is  directly  related  to  the  thickness  of  the 
equilibrium  EDL,  which  is  typically  described  as  0  (e).  Using  typical 
parameter  values  for  the  system  of  this  work,  we  calculate  a  Debye 
length  of  0.33  nm  and  e  value  of  0(10-5).  These  parameters  indicate 
a  very  small  EDL  thickness  relative  to  all  other  dimensions  of  the 
system.  We  reiterate  that  Eqs.  (1)  and  (2)  are  steady  due  to  the 
constant  boundary  conditions  and  the  liquid  phase  time  constant 
being  much  smaller  than  total  elapsed  time  during  a  capacity  test. 

Specifying  the  boundary  conditions  of  Eq.  (2)  is  a  non-trivial  task 
due  to  the  lack  of  direct  knowledge  of  the  potential  at  the  solid 
surfaces  of  each  electrode.  However  the  boundary  conditions  for 
Eq.  (1 )  are  readily  drawn  from  the  physical  consideration  that  there 
is  no  flux  of  PFg  ions  into  the  solid  phase;  that  is,  we  consider  the 
active  material  at  either  electrode  as  a  permselective  membrane 
that  only  allows  passage  of  Li+  ions.  Setting  the  PFg  flux  equal  to 
zero  gives 

^  +  cBzB ^  =0  at  x  =  0, 1  for  all  t  (4) 


+  cAzA g  =J  at  x  =  0, 1  for  all  t  (5) 

where  x  =  0  denotes  the  surface  of  the  cathode,  and  x  =  1  is  the 
coordinate  for  the  anode  surface.  Furthermore,  we  know  that  the 
flux  of  Li+  ions  is  fixed  by  the  cell  current  density  J. 

We  defer  the  discussion  of  boundary  conditions  for  Eq.  (2)  until 
later  in  this  work  because  we  show  in  the  Solution  of  the  governing 
equations  section  that  Eq.  (2)  is  not  solved  explicitly.  Instead  the 
dimensionless  form  is  used  to  justify  the  division  of  the  problem 
into  electroneutral  and  non-electroneutral  regions  and  the  result¬ 
ing  simplified  problem  is  solved  analytically. 

When  illustrating  the  definition  of  scaling  variables,  we  use 
the  (*)  superscript  to  indicate  a  dimensionless  quantity.  To 
begin  the  spatial  coordinate  x  is  scaled  as  x  =  x*/Lse P,  where 
Lsep  is  the  separator  thickness.  The  liquid  phase  concentrations  are 
scaled  based  upon  the  nominal  electrolyte  concentration 
cA  =  c*a/ca  o  and  cB  =  c*B/cBQ  where  we  note  that  cA,o  =  cB, o  =  c0 
since  the  liquid  is  macroscopically  electroneutrality.  The  classical 
definition  of  the  potential  scale  is  employed  where  </>  —  <p*{FIRT).  The 
potential  scale  is  equal  to  26  mV  at  a  temperature  of 298  K,  which  is 
the  condition  describing  all  results  in  this  work.  The  dimensionless 
current  density  J  is  defined  as  J  =  J*tfD  where  Jo  is  the  diffusion- 
limited  current  density,  Jd  =  DcoF/Lsep  which  takes  a  value  of 
approximately  630  A  m-2.  This  quantity  is  equal  to  one-fourth  the 
value  of  the  current  that  would  cause  the  electroneutral  electrolyte 
concentration  to  approach  zero  at  the  electrode  where  reduction 
occurs. 


2.2.  Electrode  solid  phase 


Fig.  l.  Modeling  domain  showing  a  separator  pore  connecting  solid  portions  Of  anode  The  diffusion  dynamics  in  the  active  material  of  each  elec- 

and  cathode  and  spatial  coordinate  definition.  Dimensions  are  not  to  scale.  trade  play  an  important  role  in  determining  the  Cell  Output 
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voltage.  Some  works  have  substituted  a  single  particle  approxi¬ 
mation  [15,16]  in  place  of  the  complexity  of  the  non-uniform 
intercalation  current  density  within  each  porous  electrode  that 
arises  from  potential  and  concentration  gradients  in  both  the 
solid  and  liquid  phases.  Our  model  is  similar  to  the  single  particle 
concept  because  a  representative  solid  diffusion  problem  is 
solved  for  each  electrode,  though  we  solve  diffusion  in  multiple 
particles  and  determine  the  local  current  for  each  particle  via  a 
voltage  constraint.  We  assume  the  intercalation  process  occurs  as 
diffusion  into  a  single  phase  in  both  graphite  and  lithium  iron 
phosphate.  This  approach  has  been  shown  to  work  well  in  a 
previous  study  [17]  and  the  exact  physics  of  the  two-phase  na¬ 
ture  of  iron  phosphate  are  still  openly  debated  [3,18],  Further¬ 
more  we  are  more  interested  in  studying  the  EDL  structure  in 
this  work  and  a  less  complex  treatment  of  the  solid  phase  is 
acceptable.  We  begin  by  describing  the  time-varying  concentra¬ 
tion  of  cations  in  the  solid  phase  of  either  electrode  in  terms  of 
the  ionic  flux.  The  convective  term  normally  associated  with 
ionic  transport  may  be  assumed  zero  if  volume  changes  in  the 
active  material  are  neglected,  and  migration  may  be  neglected  if 
significant  variations  in  potential  do  not  exist  [19],  Graphite  is  a 
good  electronic  conductor  and  therefore  the  assumption  of 
negligible  potential  variation  within  a  particle,  and  thus  negli¬ 
gible  migration,  is  applied  to  the  anode.  Though  iron  phosphate 
is  a  poor  conductor,  we  still  apply  this  assumption  to  the  cathode 
as  well  due  to  the  small  size  of  cathode  particles  leading  to 
insignificant  variation  of  potential  within  a  particle.  We  also  as¬ 
sume  a  sufficient  amount  of  electrons  are  present  to  ensure  a  net 
zero  charge  within  the  solid,  but  do  not  actively  model  their 
presence.  The  non-dimensional  solid  transport  equation  result¬ 
ing  from  the  previous  assumptions  is 


dcAJ  =  Da2cAj 

at  1  dxf 


(6) 


where  Dj  is  the  solid  phase  diffusion  coefficient  which  takes 
different  values  in  the  anode  (j  —  an)  and  cathode  (j  =  ca).  Note  that 
Xj  is  now  defined  locally  within  the  solid  phase  such  that  Xj  =  0 
refers  to  the  plane  of  symmetry  corresponding  to  the  center  of  a 
spherical  particle,  and  Xj  =  1  denotes  the  interface  between  the 
active  material  and  the  electrolyte  solvent  regardless  of  the  elec¬ 
trode.  For  the  sake  of  simplicity,  the  same  coordinate  system  is 
maintained  for  both  electrodes  (symmetry  plane  at  Xj  —  0,  solvent  at 
Xj  =  1)  even  though  in  reality  the  electrodes  would  mirror  each 
other.  Unlike  the  liquid  phase  governing  equations,  the  solid  phase 
cannot  be  assumed  to  be  in  steady  state.  Thus  a  time  scale  tj  needs 
to  be  introduced  based  on  the  characteristic  time  of  diffusion  for 
cations  within  the  solid  phase, 


where  t*  indicates  the  dimensional  time.  The  spatial  variable  is 
scaled  by  the  total  electrode  thickness 


The  concentration  scale  in  the  solid  phase  is  based  upon  the 
saturation  concentration  according  to 


.  -  A) 

1  ~  CAJ, max 


(9) 


where  Caj,max  is  the  saturation  concentration  of  the  electrode  active 
material.  The  non-dimensional  boundary  conditions  are 


dcAj 

dXj 


=  0  at  Xj  =  0, 


Scaj 

l)Xj 


=  jj  at  Xj  =  1 


(10) 


where  y j  denotes  the  volume  fraction  of  active  material  within  the 
electrode.  The  inactive  material  encompasses  conductive  additives, 
polymer  binder,  and  void  volume,  all  of  which  is  incapable  of 
allowing  lithium  intercalation.  During  the  process  of  converting  to 
non-dimensional  equations,  the  boundary  condition  becomes 
Jj  =  J*IJSj,  where  Jsj  —  DjCAj.msxF/Lj  is  analogous  to  the  diffusion- 
limited  current  density  in  the  liquid  phase  but  with  the  satura¬ 
tion  concentration  replacing  the  electroneutral  value. 

Note  that  Xj  is  defined  locally  within  the  solid  particle  such  that 
Xj  =  0  refers  to  the  center  of  the  particle  and  Xj  =  1  denotes  the 
surface.  At  the  center  of  the  particle,  there  can  be  no  flux  of  ions  due 
to  symmetry.  The  boundary  condition  at  the  interface  between  the 
particle  and  the  liquid  phase  has  a  non-zero  flux  related  to  the  local 
current  which  represents  the  intercalation  reaction  rate.  The  minus 
sign  would  be  removed  when  considering  the  cathode,  to  account  for 
the  fact  that  a  discharge  current  causes  a  flux  of  ions  into  the  active 
material.  This  information  is  summarized  in  the  schematic  of  Fig.  2. 

The  initial  conditions  of  the  solid  diffusion  problem  are  directly 
related  to  the  range  of  lithiation  incurred  by  the  electrodes,  which 
is  typically  a  design  parameter  of  Li-ion  cells  that  must  be  deter¬ 
mined  empirically.  Again  referring  to  the  conditions  for  a  capacity 
test,  we  assume  the  cell  is  fully  charged  initially  to  the  upper 
voltage  limit  of  3.6  V  which  corresponds  to  100%  SOC  as  specified  by 
the  manufacturer  [20],  A  typical  set  of  initial  non-dimensional 
degrees  of  lithiation  is 

cA,an(X;0)  =  0.80  cAa(x,  0)  =  0.025  (11) 

These  initial  conditions  are  not  known  precisely  without 
examining  experimental  data  since  the  initial  capacity  loss  during 
the  battery  formation  process  is  not  disclosed  by  the  manufacturer. 
Therefore,  the  initial  conditions  are  adjusted  empirically  until 
agreement  is  achieved  with  low  discharge  rate  experimental  data. 
It  is  noted  that  the  ionic  current  must  remain  continuous  across  the 
EDL  to  conserve  mass. 


3.  Solution  of  the  governing  equations 

This  section  details  the  solution  of  the  previously  derived  gov¬ 
erning  equations.  First  the  solution  of  the  electroneutral  liquid 
phase  is  presented  followed  by  the  electrode  solid  phase.  These 
solutions  are  utilized  by  the  polarographic  relationships  in  the 
Helmholtz  and  Gouy— Chapman  models  to  predict  the  time-varying 
cell  voltage. 

3.1.  Solutions  in  the  electrode  solid  phase 

An  analytical  solution  is  obtainable  for  a  steady  flux.  Normally  in 
non-homogeneous  problems,  a  steady  solution  would  be  defined 


Xj=  0  Xj=Lj 


Fig.  2.  Summary  of  the  boundary  conditions  for  solid  diffusion  within  electrode  active 
material,  presented  in  general  form  for  either  the  anode  or  the  cathode. 
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that  satisfies  the  corresponding  boundary  conditions.  However, 
with  the  prescribed  boundary  conditions  of  this  problem,  there  is 
no  steady  solution  because  the  salt  concentration  within  the 
modeling  domain  continually  increases  in  time.  Instead,  the  solu¬ 
tion  is  broken  into  three  parts  [21  ]  corresponding  to  a  steady  term 
that  defines  the  shape  of  the  long-time  response,  but  not  its 
magnitude;  a  transient  term  that  defines  the  difference  between 
the  long-time  response  and  the  actual  solution  at  time  t;  and  the 
time-varying  mean  concentration. 

First  the  solution  for  the  mean  concentration  is  computed  by 
integrating  the  non-dimensional  form  of  Eq.  (6)  and  applying  the 
boundary  conditions  of  Eq.  (10).  Then,  the  steady  profile  term  is 
computed  based  on  the  boundary  conditions  and  time-varying 
mean  concentration.  Finally,  the  transient  term  is  solved  by  sepa¬ 
ration  of  variables.  Doing  so  yields 


cAj(x,t)  =  CAj(X,  0)  +Jj  I  tj  H 


3x?-l 


-^2  E  (-^Te  nVtjC0S 


The  first  one  hundred  terms  of  the  infinite  summation  are  used 
when  computing  the  solution.  If  the  external  current  density  is 
fixed  then  the  solution  for  the  bulk,  steady  state,  liquid  concen¬ 
tration  and  potential,  as  well  as  the  time-varying  solid  concentra¬ 
tion  in  each  electrode,  is  given  for  all  time.  For  t  much  greater  than 
unity,  the  solution  becomes  linearly  increasing  with  respect  to  t  and 
parabolic  in  the  spatial  coordinate  x. 

3.2.  Solutions  in  the  electroneutral  liquid 

At  this  point,  we  wish  to  derive  an  analytical  solution  for  the 
concentration  and  potential  within  the  separator  defined  as 
0(1)  —  0(0).  Eqs.  (1)  and  (2)  with  corresponding  boundary  condi¬ 
tions  comprise  a  system  of  three  coupled  ordinary  differential 
equations  when  considering  the  ionic  concentrations  indepen¬ 
dently.  The  coupling  is  stiff  due  to  the  introduction  of  the  e 
parameter.  In  order  to  make  the  system  mathematically  tractable, 
we  divide  our  analysis  to  deal  with  the  electroneutral  region 
(commonly  referred  to  as  the  ’bulk’)  and  the  EDL  separately.  To 
develop  an  equation  for  the  bulk  solution  we  first  note  that  in  the 
bulk  the  ionic  concentrations  are  equal  (ca  (x)  =  Cb(x)  =  c(x)) 
everywhere  to  establish  electroneutrality.  Then  we  sum  the  equa¬ 
tions  for  each  individual  ionic  concentration,  in  non-dimensional 
form,  and  integrate  with  J  constant  to  obtain 


J  =  2 


0C 


To  obtain  an  equation  for  the  electric  potential  the  equation 
describing  the  ionic  concentration  of  the  anion  is  subtracted  from 
that  of  the  cation,  and  the  result  is  integrated  to  yield 


Since  an  expression  for  the  bulk  concentration  has  already  been 
derived  as  Eq.  (15),  it  may  be  substituted  into  Eq.  (16)  and  inte¬ 
gration  performed  again  to  obtain 


0(X) 


-HH)) 


(17) 


(12) 


The  solution  presented  in  Eq.  (17)  is  logarithmic  because  when 
current  is  passed  through  an  electrolyte,  a  small  volume  charge 
density  of  0(e2)  is  generated.  This  is  not  expected  to  satisfy  Eq.  (2) 
because  the  solution  phase  is  only  approximately  electrically 
neutral. 

3.3.  Solutions  in  the  electrical  double  layer 

The  solution  for  potential  in  the  EDL  is  found  by  solving  Pois¬ 
son’s  equation  subject  to  the  Boltzmann  distribution  of  ionic  con¬ 
centration,  with  the  additional  boundary  condition  imposed  by  the 
total  potential  calculation  of  the  reaction  kinetics. 


=  sinh(0EDL) 


d20 

d y2 

To  solve  Eq.  (18),  multiply  by  2(d0/dy)  and  integrate  twice  to 
obtain  the  inner  potential 


0EDLj(y)  =  4  tan/i  1  ftanh^e 


(19) 


where  y  —  x/e,  Zj  represents  the  diffuse  layer  potential,  also 
referred  to  as  the  zeta  potential,  and  again  all  potentials  are  scaled 
by  RT/F.  The  boundary  conditions  for  the  potential  of  Eq.  (19)  are 
0EDL,ca  —  Cca  at  x  —  0  and  0EDL,an  =  Can  at  x  =  1.  The  additional 
required  matching  conditions  are 

ylim  0EDL,ca  O')  =  Hin  0(X)  at  X  =  0,  ^lim  0EDL,an  (y) 

=  lint  0(x)  at  x  =  1,  (20) 

The  concentration  solution  in  the  EDL  relies  on  the  Boltzmann 
distribution  for  ionic  concentration  which  is  valid  for  equilibrium 
conditions 


(13)  q(y)  =  tmJr****M 


We  will  solve  Eq.  (13)  directly  via  integration.  The  constant  of 
integration  must  be  determined  by  applying  conservation  of  mass 
for  the  anion 

I  cB(x)  dx  =  1  (14) 

0 

This  conservation  of  mass  condition  comes  from  the  constraint 
that  there  is  no  flux  of  PFg  ions  into  the  solid  portion  of  either 
electrode,  which  implies  that  the  total  number  of  anions  within  the 
liquid  must  be  constant.  Integrating  Eq.  (13)  and  applying  Eq.  (14) 


where  c^j  is  the  concentration  at  the  interface  between  the  EDL 
and  the  electroneutral  liquid  as  x  tends  to  zero  or  one.  Eqs.  (21 )  and 
(19)  are  collectively  referred  to  as  the  Poisson-Boltzmann  distri¬ 
bution  for  the  classical  EDL. 

3.4.  Polarographic  relationships 

The  entire  liquid  phase  solution,  encompassing  solutions  in  the 
electroneutral  portion  as  well  as  the  EDL,  may  be  viewed  in  terms  of 
a  polarographic  relationship  for  the  case  of  a  steady  external  cur¬ 
rent.  The  term  polarographic  refers  to  the  steady  state  relationship 
between  current  and  voltage  which  defines  the  amount  of  over¬ 
potential  required  to  sustain  a  given  current  through  the  liquid 
phase,  and  this  overpotential  manifests  itself  as  a  deviation  from 
the  cell  open  circuit  voltage. 
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Before  any  modeling  of  reaction  kinetics  can  begin  the  exact 
chemical  reactions  being  studied  must  be  specified.  For  this  dis¬ 
cussion  we  consider  the  processes  that  occur  during  discharge, 
though  the  processes  that  occur  during  charge  are  not  difficult  to 
determine  simply  by  reversing  the  direction  of  current  flow  in  the 
cell  and  following  the  order  of  events  backwards.  In  the  solid  phase 
of  the  anode  lithium  exists  in  its  atomic  form.  When  a  discharge 
current  is  present,  lithium  de-intercalates  from  the  host  matrix  and 
enters  the  liquid  phase,  losing  an  electron  in  a  process  known  as 
oxidation  and  becoming  a  lithium  ion  with  valence  equal  to  one. 
This  process  creates  an  open  site  in  the  host  matrix  which  the 
lithium  atom  once  occupied.  For  the  chemistry  of  interest,  the  re¬ 
action  is  given  by 

LiaC6  <->Li+  +  e~  +  6C  (22) 

where  atomic  lithium  in  the  anode  solid  phase  (c/iian)  is  represented 
by  the  left-hand  side  term  LiaC6,  the  lithium  ions  in  the  liquid  phase 
(ca)  are  represented  by  Li+,  and  the  open  site  in  the  anode  solid 
phase  created  by  de-intercalation  (c4ian,max  -  Ca,an)  is  given  by  6C.  In 
the  cathode  during  discharge,  the  lithium  ions  from  the  liquid 
phase  intercalate  into  the  solid  phase.  In  doing  so,  they  gain  an 
electron  in  a  process  referred  to  as  reduction  and  thus  again 
become  atomic  lithium  with  no  charge.  Additionally,  an  open  site 
that  was  once  available  is  now  occupied  by  a  lithium  atom.  The 
reaction  is  defined  by 

LibFeP04<->Li+  +  e~  +  FeP04  (23) 

where  the  lithium  ions  in  the  liquid  phase  are  again  represented  by 
Li+,  the  open  site  within  the  cathode  solid  phase  (c4  ca,max  -  c4iCa)  is 
given  by  FeP04,  and  the  atomic  lithium  in  the  cathode  solid  phase 
(c4,Ca)  is  defined  by  LibFeP04.  These  reactions  govern  the  charge/ 
discharge  process  with  the  right  arrow  indicating  oxidation  and  the 
left  arrow  indicating  reduction.  We  introduce  the  convention  that 
oxidation  corresponds  to  a  positive  reaction  current  and  reduction 
to  a  negative  reaction  current  in  both  electrodes. 

Le  Chatelier’s  principle  states  that  when  a  chemical  system  at 
equilibrium  experiences  a  change  in  operating  conditions  such  as 
reactant  concentrations  or  temperature,  the  equilibrium  will  shift 
to  counteract  the  change  and  a  new  equilibrium  will  be  established 
[22],  We  can  apply  this  principle  to  Eqs.  (22)  and  (23)  to  derive  an 
appropriate  form  of  the  Butler-Volmer  kinetic  expressions 
defining  the  relationship  between  overpotential,  temperature, 
concentration,  and  the  intercalation  rate.  An  increase  in  LiaC6  will 
cause  the  rate  of  the  oxidation  reaction  to  increase,  whereas  an 
increase  in  Li+  or  6C  will  cause  the  rate  of  the  reduction  reaction  to 
increase.  For  Eq.  (23),  an  increase  in  Li+  or  FeP04  will  increase  the 
rate  of  the  oxidation  reaction  while  an  increase  in  LibFeP04  will 
increase  the  rate  of  the  reduction  reaction. 

A  brief  discussion  of  why  the  quantities  c4jimax  —  c4j  are  intro¬ 
duced  is  warranted.  Considering  the  anode,  the  symbol  6C  repre¬ 
sents  an  open  site  in  the  solid.  The  number  of  open  sites  is  equal  to 
the  total  number  of  sites  per  unit  volume,  c4,an,max,  minus  the 
number  of  occupied  sites  per  unit  volume,  c4  an.  Therefore  the 
symbol  6C  is  equivalent  to  the  quantity  c4  an  max  -  cAan,  and  a 
similar  discussion  applies  in  the  cathode  for  the  relationship  be¬ 
tween  c4iCa,max  -  Ca, ca  and  the  symbol  FeP04.  The  rate  equations 
resulting  from  the  reaction  definition  and  analysis  invoking  Le 
Chatelier’s  principle  [22]  are 

Jj  =  KA/&***  -  ~  (M) 

where  each  k*  is  a  dimensional  reaction  rate  constant  which 
relates  the  conditions  of  the  reaction  to  the  net  reaction  rate, 


A0*j  is  the  Stern  layer  potential  difference  which  will  be  dis¬ 
cussed  in  greater  detail  in  the  following  sections,  and 
^  oo  ca  =  ^(0)  for  the  cathode  and  cA  ^  an  =  c*A(\ )  for  the  anode 
is  the  concentration  of  Li+  ions  at  the  interface  between  the  EDL 
and  the  electroneutral  liquid.  The  parameter  a  is  referred  to  as  a 
symmetry  factor  or  transfer  coefficient.  It  defines  the  fraction  of 
applied  potential  that  contributes  to  the  oxidation  (1  -  a)  and 
reduction  (a)  reaction  rates  and  is  typically  taken  to  be  a  —  0.5 
due  to  the  symmetry  exhibited  by  most  experimental  current- 
voltage  relationships  with  respect  to  charge  and  discharge.  These 
rate  equations  will  be  used  to  derive  an  analytical  current- 
voltage  relationship  for  a  steady  current.  They  are  similar  to  the 
relationships  presented  by  Bazant  ],  but  with  the  addition  of 
the  dependence  on  filled  and  open  sites  in  the  active  material 
matrix. 

As  stated  previously,  the  Stem  layer  voltage  is  only  part  of  the 
entire  potential  drop  within  the  EDL  in  the  general  case.  In  fact  the 
entire  potential  difference  of  the  EDL,  A0eDLj,  is  equal  to  the  sum  of 
the  Stern  layer  and  diffuse  layer  potentials  such  that 
A0edlj  =  A <pSj  +  ‘Cj.  Visually,  we  represent  the  assignment  of  the 
total  EDL  potential  difference  to  the  Stern  layer  and  the  diffuse  layer 
in  Fig.  3. 

A  positive  current  with  J  >  0,  directionally  from  the  anode  to  the 
cathode,  requires  a  positive  potential  difference  across  the  anode 
EDL  (A$EDL,an  =  $ an  -  $(1),  A$EDL,an  >  0)  and  a  negative  potential 
difference  across  the  cathode  EDL  (Afeuca  =  4>a  -  $(0), 
A0EDL,ca  <  0).  Imposing  J  >  0  also  corresponds  to  discharging  the 
cell  in  the  sign  convention  of  this  work,  which  requires  Jan  >  0  and 
Jca  <  0  to  match  the  oxidation/reduction  sign  convention.  From  this 
point  forward,  we  omit  subscripts  unless  absolutely  necessary  and 
note  that  for  a  positive  (discharge)  current  Jan  =  J  and  Jca  =  -J  due  to 
the  sign  convention  associated  with  oxidation  and  reduction 
currents. 

A  non-dimensional  form  of  Eq.  (24)  is  obtained  by  introducing 
the  current  scale,  Jsj,  the  potential  scale,  and  the  concentration 
scales  of  the  solid  and  liquid.  Then,  the  rate  constants  are  re-defined 
as  koj  =  fco/maxj/Jsj  and  krj  =  k*jC0cmaxj/JSj,  to  avoid  carrying 
several  constants  continuously.  The  non-dimensional  rate  equation 
that  results  is 


Jj  =  koJCAA-“^  -  krjcA,„j(  1  -  cAj)e-^  (25) 

The  pre-multipliers  of  the  exponential  terms  are  directly  related 
to  the  state  of  charge  of  each  electrode.  As  stated  in  the  introduction 


x  =  0  x  x  =  1 


Fig.  3.  Definition  of  potential  differences  in  the  electrical  double  layer  at  anode  and 
cathode.  The  potential  profile  in  the  bulk  is  omitted. 
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this  provides  a  physical  mechanism  for  altering  the  values  of  the 
rate  constants  during  the  discharge  process,  or  in  this  instance  the 
product  of  the  rate  constant  and  the  reactants. 


where  y  =  x/e  and  the  common  part  is  equal  to  the  bulk  value  at 
x  =  0  and  x  =  1.  Following  a  similar  methodology,  the  complete 
solution  for  concentration  is 


3.4.1  Helmholtz  model 

We  note  that  in  general  the  voltage  difference  across  the  sepa¬ 
rator  (referred  to  interchangeably  as  the  liquid  voltage,  14)  is  given 

by 


q(x)  =  1  +{(*  ~l)  +  (l  -  l) 

+  (1+i)(^)-1) 


(30) 


14(t)  =  A 0gDL  an  +  2  tanh1  (J)  +  d0EDLca  (26) 

where  the  inverse  hyperbolic  tangent  term  arises  from  the  poten¬ 
tial  difference  in  the  electroneutral  liquid,  and  the  other  two  terms 
are  associated  with  the  potential  difference  across  the  electrical 
double  layer.  Upon  solving  for  the  Stern  layer  potential  of  Eq.  (24), 
non-dimensionalization,  and  substitution  into  Eq.  (26),  the  current 
voltage  relationship  is  derived  as 


where  <p(x)  is  defined  in  Eq.  (29).  Noting  that  fe)L,ca  —  -?ca  due  to 
the  convention  of  Fig.  3  where  the  zeta  potential  is  defined  relative 
to  the  bulk  liquid  potential,  and  substituting  the  potential  from  Eq. 
(24),  the  current  voltage  relationship  of  Eq.  (26)  becomes 


14(f)  =  4  tanh-1  (J)  +  In 


\an(l  ~  CAan) 
<  ^anCAan -J 


kp,caCA,ca  +J  \ 

kr,ca(l-CAca); 


(31) 


v(j+  sjj2  +4k0,anC^ankrianCA(l)(l-C/,,an)  ^  ^  (27) 

X  \-J  +  V/i2+4ko,caCAcakr,caCA(0)(l  -CAca)  J  J 


The  first  term,  2  tan  h  1  (/),  defines  the  voltage  in  the  electro¬ 
neutral  liquid.  The  natural  logarithm  term  represents  the  voltage  of 
the  EDLs  at  the  anode  and  cathode. 

The  current-voltage  relationship  must  be  modified  to  account 
for  the  use  of  an  equilibrium  potential  during  model  develop¬ 
ment.  Typical  open-circuit  voltages  for  the  anode  and  cathode 
are  empirically  obtained  using  half-cell  measurements,  where 
the  working  electrode  is  assessed  relative  to  a  lithium  metal 
counter/reference  electrode.  The  zero  current  EDL  potential  is 
inherently  included  in  the  measured  data,  so  the  EDL  voltage 
value  at  zero  current  is  subtracted  from  the  measured  open- 
circuit  voltage.  The  overall  process  is  similar  for  either  the 
Helmholtz  or  the  Gouy— Chapman  models  of  the  following  sec¬ 
tion,  with  the  only  difference  being  a  change  in  the  form  of  the 
polarographic  expression.  Mathematically,  the  preceding  dis¬ 
cussion  is  accomplished  as 


The  electroneutral  portion  of  the  EDL  voltage,  2  tan  ft-1  (J),  re¬ 
mains  as  discussed  previously.  However  an  additional  2  tan  h-1  (/) 
and  natural  logarithm  term  are  associated  with  the  EDL  voltage  in 
the  Gouy-Chapman  model.  Again  Eq.  (28)  must  be  applied  to  ac¬ 
count  for  only  the  non-equilibrium  portion  of  the  EDL  in  the  cell 
voltage.  This  is  different  from  the  expression  derived  using  the 
assumption  that  the  entire  potential  difference  of  the  EDL  occurs  in 
the  Stern  layer.  We  will  further  explore  the  similarities  and  differ¬ 
ences  between  Eqs.  (27)  and  (31)  in  the  Results  and  discussion 
section,  including  comparing  models  using  each  hypothesis  to 
experimental  data.  The  aim  of  this  comparison  is  to  understand  if 
one  model  more  accurately  describes  the  nano-scale  physics  of  a 
commercially  available  Li-ion  cell.  As  discussed  in  Refs.  [10],  there  is 
an  infinite  number  of  possibilities  for  division  of  the  entire  EDL 
potential  across  the  Stern  layer  and  the  diffuse  layer  determined  by 
the  capacitance  of  each  layer.  However  we  note  in  the  Results  and 
discussion  that  the  results  are  ultimately  quite  similar  when 
comparing  the  two  extremes  discussed  here,  so  further  investiga¬ 
tion  of  the  intermediate  cases  is  outside  the  scope  of  this  paper.  An 
interesting  result  of  the  preceding  sections  is  that  the  polarization 
required  in  the  EDL  is  determined  in  closed  form  by  the  battery 
SOC,  the  current,  and  the  solution  of  the  fully  electroneutral  portion 
of  the  liquid  in  both  cases. 


4.  Results  and  discussion 


VL(t)  =  VL(t)  -  (28) 

where  14  (t)  on  the  right-hand  side  is  defined  in  Eq.  (27)  and  V4(t) 
on  the  left-hand  side  accounts  for  only  the  non-equilibrium  portion 
of  the  EDL  potential. 

3.4.2.  Gouy— Chapman  model 

Before  presenting  the  polarographic  relationship,  the  complete 
solutions  for  potential  and  concentration  are  given.  The  complete 
solution  for  potential  is  equal  to  the  electroneutral  solution  given 
by  Eq.  (17),  plus  the  inner  solution  defined  by  Eq.  (19),  minus  their 
common  part 

^--n((l(;)j>1++;))+4--,H(%)^ 

+  4  tanh-1  ^tanh  e-(’-y)^ 

(29) 


In  this  section,  we  present  the  results  obtained  from  the  solution 
of  the  governing  equations  for  a  galvanostatic  capacity  test.  We  also 
discuss  the  applicability  of  the  two  descriptions  of  the  EDL  in  view 
of  system-level  considerations  and  experimental  data,  and  finally 
introduce  empirical  modifications  to  the  theoretical  current- 
voltage  relationship. 

4.1.  Definition  of  cell  terminal  voltage 

The  output  voltage  is  obtained  by  moving  from  cathode^  to 
anode  and  summing  the  voltages  encountered,  recalling  that  VL(t) 
contains  the  effects  of  both  the  electroneutral  liquid  as  well  as  the 
EDL 

V*(t)  =  U*a(cca(L,  t))  -  Uan(can(L,  t))  -  VL(t)  -  RJ*(t)  (32) 

An  ohmic  resistance  Rc  accounts  for  the  initial  ohmic  polariza¬ 
tion  that  is  not  due  to  the  bulk  liquid.  Essentially  it  is  the  sum  of 
ohmic  losses  from  the  solid  phase  of  the  electrode,  the  solid- 
electrolyte  interphase  (SEI)  layer,  and  poor  contact  between  iron 
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phosphate  particles  and  the  conductive  matrix  [3,17],  The  sum  of 
voltage  losses  is  treated  as  a  deviation  away  from  the  open-circuit 
voltage  between  electrodes  given  by  the  difference  of  U*3  and  U3n. 
The  difference  in  the  open-circuit  voltages  is  taken  as  the  ther¬ 
modynamic  maximum  voltage  that  can  be  achieved  under  ideal 
conditions  during  discharge  [23],  though  the  concentration  over¬ 
potential  in  the  solid  phase  is  embedded  in  these  functions  through 
their  dependence  on  the  surface  concentration  rather  than  the 
mean  solid  concentration.  These  functions  are  defined  in  the 
Appendix,  and  they  must  be  obtained  empirically  due  to  the  de¬ 
viations  away  from  the  ideal  case  represented  by  the  Nemst 
equation  [24],  Essentially  Eq.  (32)  states  that  the  thermodynamic 
equilibrium  voltage  is  modified  to  account  for  losses  arising  from 
concentration  polarization  in  the  solid  and  liquid,  charge  transfer, 
and  ohmic  resistances. 

The  parameter  values  used  to  produce  the  simulation  results  of 
the  following  section  are  documented  in  Table  1.  The  source  of  each 
parameter  is  defined  by  the  superscript  next  to  the  parameter 
definition.  For  measured  parameters,  AFM  and  SEM  images  from 
the  Center  for  Automotive  Research  and  the  Nanoprobe  Laboratory 
for  Bio-  and  Nanotechnology  and  Biomimetics  were  examined  us¬ 
ing  image  analysis  software.  Active  material  volume  fractions  were 
tuned  based  on  the  measured  electrode  volumes  to  match  the 
composition  range  of  the  electrodes  given  by  Ref.  [17],  Since  the 
anode  is  the  capacity-limiting  electrode  during  discharge  in  this 
instance,  its  diffusion  coefficient  was  estimated  using  the  capacity 
difference  between  capacity  tests  of  varying  current  magnitude. 
Specifically,  the  value  was  adjusted  to  account  for  the  difference  in 
capacity  between  the  C/3  and  C/1.2  tests  reported  in  the  Experi¬ 
mental  section.  Finally  the  rate  constants  are  adjusted  using  the 
following  process.  First  the  experimentally  measured  voltage  for 
the  C/3  rate  is  subtracted  from  the  cell  open-circuit  voltage  defined 
by  the  first  two  terms  of  Eq.  (32).  This  gives  the  amount  of  over¬ 
potential  experienced  by  the  cell  over  the  course  of  a  discharge 
cycle.  Then  the  rate  constant  values  and  contact  resistance  are 


Table  1 

Summary  of  model  parameters  used  for  simulation  results. 


Parameter 

Definition 

Value 

Tan 

Anode  active  material  volume  fraction11 

0.44 

Tea 

Cathode  active  material  volume  fractionb 

0.36 

Liquid  permittivity' 

46  F  m-1 

A 

Current  collector  area3 

0.19  m2 

Cajnaxan 

Saturation  concentration  of  Li+ 

30,500  mol  m~3 

in  anode  active  material' 

Saturation  concentration  of  Li+ 

16,300  mol  m-3 

in  cathode  active  material' 

Ca„(Xan,0) 

Initial  degree  of  anode  intercalation' 

0.80 

Cca(Xca,0) 

Initial  degree  of  cathode  intercalation' 

0.025 

Electroneutral  concentration' 

1000  mol  m-3 

D 

Effective  liquid  phase  diffusion  coefficient' 

1.7  x  10-10  m2  s"1 

Dan 

Anode  solid  phase  diffusion  coefficient11 

4.0  x  10“14  m2  s-1 

Dca 

Cathode  solid  phase  diffusion  coefficient' 

8.0  x  10-18  m2  s-1 

tan 

Anode  effective  diffusion  length*1 

3.4  (im 

La 

Cathode  effective  diffusion  length*1 

31  nm 

Lsep 

Separator  thickness' 

25  (im 

Jcoan 

Anode  oxidation  rate  constant*1 

2.5  x  10~3  A  m 

mol-1  s-1 

fco,ca 

Cathode  oxidation  rate  constant*1 

4.0  x  10-7  Am 

mol-1  s-1 

L.a„ 

Anode  reduction  rate  constant*1 

2.5  x  10“3  A  m4 

kr.ea 

Cathode  reduction  rate  constant*1 

4.0  x  10-7  A  m4 

mol-2  s-1 

Re 

Contact  resistance*1 

26  x  10-3  Q  m2 

T 

Temperature 

298  K 

Sources  are:  a)  measured,  b)  estimated  from  available  experimental  data,  c)  from 
literature. 


adjusted  so  that  the  last  two  terms  of  Eq.  (32)  give  an  adequate  fit  to 
the  voltage  losses.  Since  there  are  four  values  of  the  rate  constant, 
two  assumptions  are  made  to  simplify  the  fitting  process.  First,  the 
forward  and  backward  rate  constants  are  assumed  to  be  equal, 
since  there  is  no  evidence  to  the  contrary  in  typical  intercalation 
batteries.  Second,  the  anode  rate  constant  values  are  set  high 
enough  to  give  minimal  contribution  to  the  overpotential,  since  this 
is  generally  the  case  for  graphite  electrodes  operating  near  room 
temperature.  This  leaves  only  the  cathode  rate  constant  and  contact 
resistance  to  fit  the  observed  overpotential.  The  rate  constant 
values  are  identified  only  for  the  lowest  current  and  then  used  to 
predict  the  voltage  at  higher  currents.  To  add  another  set  of  data  to 
test  the  predictive  nature  of  the  model,  3C  capacity  test  data  from 
Ref.  [25]  is  plotted  as  well. 

4.2.  Analytical  predictions 

The  steady  state  profiles  for  concentration  using  each  of  the  EDL 
hypotheses  are  plotted  in  Fig.  4.  The  profiles  computed  using  the 
Gouy— Chapman  model  will  be  largely  the  same  as  those  for  the 
Helmholtz  model,  except  for  a  small  region  of  0(10-5)  near  either 
solid  boundary.  Since  the  liquid  phase  transport  reaches  steady- 
state  quickly  compared  to  the  total  discharge  time,  the  concentra¬ 
tion  in  the  liquid  is  a  function  of  only  space  and  current  density.  As 
the  limiting  discharge  current  (J  =  4)  is  approached  the  concen¬ 
tration  at  the  cathode  surface  goes  to  zero.  During  a  charge  con¬ 
dition,  the  steady-state  concentration  near  the  anode  would 
approach  zero  for  currents  close  to  the  limiting  value. 

The  solution  for  potential  is  plotted  in  Fig.  5,  showing  it  is  log¬ 
arithmic  with  respect  to  the  spatial  coordinate.  The  total  potential 
difference  is  non-linear  with  respect  to  the  current  and  this  non¬ 
linearity  influences  the  concentrations  at  the  Stern  plane  as  well. 
The  potential  difference  between  the  electrodes  approaches  zero  as 
the  current  approaches  its  limiting  value.  From  cell  level  mea¬ 
surements  one  could  expect  a  large  magnitude  of  overpotential  for 
currents  near  the  limit,  though  in  reality  this  limit  is  not  easily 
approached  before  other  limitations  such  as  ohmic  overpotential  or 
solid  diffusion  become  prohibitive. 

Next  the  solid  phase  diffusion  solution  is  examined.  The  spatial 
distribution  of  Li+  concentration  within  the  solid  phase  of  each 
electrode  is  plotted  at  various  times  during  a  discharge  at  4.8C  rate 
in  Fig.  6.  This  current  is  chosen  because  it  is  large  enough  to 
highlight  the  spatial  gradients  that  occur,  which  tend  to  be  less 
apparent  at  lower  rates.  Since  the  cathode  actually  exhibits  a  phase 
change,  solving  the  diffusion  equation  is  an  empirical  method  of 
accounting  for  rate  effects. 

4.3.  Polarographic  relationship  analysis 

Fig.  7a  gives  a  comparison  of  the  polarographic  expression 
relating  current  and  steady  state  voltage  in  the  liquid  for  the 
Helmholtz  and  Gouy-Chapman  models.  Recall  that  this  is  a  sum  of 
the  potential  difference  within  the  electroneutral  liquid  and  the 
overpotential  of  the  EDL  To  produce  Fig.  7a,  the  solid  concentration 
values  are  set  to  their  corresponding  values  for  each  SOC  of  interest. 
Then  the  current  is  swept  from  zero  to  the  limiting  value  while 
holding  all  other  variables  constant  to  give  a  polarographic  curve 
for  the  liquid  between  electrodes.  The  plotted  curves  resulting  from 
each  limiting  description  of  the  EDL  are  fairly  similar  in  several 
respects.  First,  both  curves  tend  to  shift  upward  or  downward  as  the 
cell  SOC  is  varied.  However  the  Helmholtz  model  exhibits  sym¬ 
metry  about  50%  SOC  in  its  polarographic  curves  whereas  the 
Gouy— Chapman  model  does  not.  The  polarographic  curves  of  the 
Gouy— Chapman  model  also  tend  to  shift  by  a  larger  amount  as  the 
cell  state  of  charge  varies. 
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Another  interesting  point  of  Fig.  7a  is  the  voltage  response  at 
low  current  density.  Due  to  the  high  value  of  the  diffusion  limited 
current  density  in  this  case,  we  would  not  expect  to  see  operation  in 
the  region  [/|  >  0.5.  Furthermore  for  current  rates  of  less  than  or 


equal  to  1C,  we  have  [/|  <  0.02,  and  this  corresponds  to  the  region 
displayed  in  the  inset  of  Fig.  7a.  In  the  context  of  lithium  ion  battery 
literature,  the  double  layers  are  assumed  to  be  thin  and  the  Butler— 
Volmer  kinetic  law  is  applied  based  on  concentration  and  potential 


(a) 

Fig.  6.  (a)  Concentration  of  Li+  within  (a)  cathode  and  (b)  anode  given  by  Eq.  (12)  for  a  full 
positioned  at  Xj  =  0  and  interface  between  active  material  and  solvent  is  at  Xj  =  1.  Time  : 


(b) 

I  discharge  cycle  with  J  =  8.7  x  10  2  and  parameter  values  from  Table  1.  Symmetry  plane  is 
scale  shown  is  based  on  the  anode. 
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values  in  the  bulk.  Reflecting  on  our  solutions,  this  essentially  im¬ 
plies  that  battery  literature  employs  the  Helmholtz  model  for  the 
EDL  and  neglects  the  possibility  of  diffuse  charge  caused  by  a  non- 
negligible  zeta  potential. 

To  continue  the  comparison  between  the  two  models,  the 
overpotential  due  to  charge  transfer  in  the  Helmholtz  model  is 
plotted  in  Fig.  7b  along  with  the  overpotential  according  to  the 
Gouy-Chapman  model.  The  comparison  is  made  at  the  C/3,  C/1.2, 
and  4.8C  rates  which  corresponds  to  J  =  5.6  x  10-3, 1.5  x  10-2,  and 
8.7  x  10-2  respectively.  This  differs  from  the  results  of  Fig.  7a 
because  here  we  are  fixing  a  current  and  computing  the  resulting 
time-varying  response  of  concentration  and  potential.  We  note  that 
the  curves  show  qualitative  similarity  but  the  Gouy-Chapman 
model  exhibits  significantly  greater  variation  in  voltage  over  the 
course  of  discharge.  Additionally  the  Gouy-Chapman  model  at  the 
highest  current  level  shows  the  beginning  of  an  increase  in  over¬ 
potential  near  the  end  of  discharge  that  could  be  associated  with  a 
voltage  knee  if  it  were  larger  in  magnitude. 

Fig.  8a  shows  the  polarographic  curves  that  result  from  several 
different  choices  of  the  kinetic  rate  constants.  The  selection  of  the 
rate  constants  is  an  uncertain  process  that  requires  comparisons 
with  experimental  data.  The  EDL  voltage  varies  nonlinearly  as  the 


rate  constants  are  changed.  Lower  rate  constants  lead  to  higher  EDL 
voltage,  and  higher  rate  constants  lead  to  lower  EDL  voltage.  With 
higher  rate  constants,  less  overpotential  is  required  to  pass  a  given 
amount  of  current. 

Ref.  [17]  has  noted  difficulty  in  choosing  a  rate  constant  for  the 
LiFePCU  system.  In  particular,  they  found  that  rate  constants  that  fit 
low  rate  experimental  data  overpredicted  the  overpotential  at 
higher  rates.  To  mitigate  this  issue,  they  defined  the  rate  constant  as 
an  empirical  function  of  the  total  cell  current  density.  Fig.  8b  gives 
an  illustration  of  why  this  approach  is  not  repeated  here.  The 
overpotential  predicted  by  the  polarographic  expressions  of  Eqs. 
(27)  and  (31)  reasonably  approximates  the  literature  prediction 
which  uses  a  variable  rate  constant  that  is  a  function  of  the  current 
density.  The  comparison  cannot  be  extended  to  higher  rates, 
because  the  polynomial  fit  used  in  Ref.  [17]  has  an  inflection  point 
beyond  the  presented  current  densities  and  eventually  becomes 
negative. 

5.  Experimental 

Data  obtained  from  a  cylindrical  graphite/iron  phosphate  cell 
with  can  dimensions  of  26  mm  diameter  by  65  mm  height  (26,650) 


Fig.  8.  (a)  The  EDL  voltage  as  a  function  of  current  for  selected  cathode  rate  constant  values.  The  anode  rate  constants  are  fixed  to  2.5  x  10“3.  Units  of  the  oxidation  rate  constants 
are  A  m  mol  1  s  ',  and  units  of  the  reduction  rate  constants  are  A  m4  mol-2  s  (b)  Comparison  of  the  literature  values  for  kinetic  overpotential  using  a  variable  rate  constant  with 
the  output  of  Eqs.  (27)  and  (31).  using  parameter  values  of  Table  1. 
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are  used  for  comparative  purposes  and  to  examine  the  validity  of 
the  double-layer  assumptions.  Galvanostatic  capacity  tests  are 
performed  at  current  rates  of  C/3,  C/1.2,  and  4.8C.  Data  is  collected 
at  10  Hz  and  current  is  commanded  using  a  PLA800-60-300  power 
load  and  supply  from  American  Reliance,  Inc.  Temperature  is  fixed 
at  25  °C  using  an  AC-027  Peltier  junction  from  TE  Technology,  Inc. 
along  with  sufficient  insulation  to  prevent  excessive  power  draw  by 
the  Peltier  junction.  The  lower  voltage  limit  is  2.5  V  while  the  upper 
limit  is  3.6  V,  as  defined  by  the  manufacturer  [20],  During  charging, 
constant  voltage  is  enforced  at  the  upper  voltage  limit  for  30  min. 

5.1.  Comparison  of  theoretical  model  structure  with  experimental 
data 

The  model  output  using  the  Helmholtz  hypothesis  is  compared 
with  experimental  data  in  Fig.  9a  for  multiple  current  rates.  We  see 
there  is  generally  good  agreement  between  the  steady-state  model 
using  the  Helmholtz  model  of  the  EDL  and  the  experimental  data 
gathered  at  low  current.  The  model  results  diverge  from  the 
experiment  for  the  highest  discharge  rates.  This  is  to  be  expected 
however,  because  several  factors  are  neglected  in  developing  this 
simplified  model  that  would  have  a  growing  impact  as  the  current 
rate  increases.  These  factors  include  the  neglect  of  the  non-uniform 
reaction  rate  with  respect  to  the  thickness  of  the  electrode,  neglect 
of  thermal  effects,  and  neglect  of  the  resistive  reactant  nature  of 
lithium  iron  phosphate  [4,17],  Given  that  the  model  overpredicts 
the  cell  voltage  compared  to  the  experimental  data  at  the  highest 
current,  the  most  likely  factor  to  improve  the  model  predictions  is 
inclusion  of  the  time-varying  resistance  of  the  cathode  that  results 
from  the  resistive  reactant  effects. 

The  model  output  using  the  Gouy-Chapman  hypothesis  is 
compared  with  experimental  data  in  Fig.  9b  for  multiple  current 
rates.  Recall  that  the  Gouy— Chapman  hypothesis  allowing  for 
mobile  charge  in  the  diffuse  layer  gives  a  wider  variation  in  the 
kinetic  overpotential  over  the  course  of  the  discharge.  The  wide 
variation  in  the  overpotential  magnitude  leads  to  the  non-physical 
phenomena  of  the  voltage  actually  increasing  during  the  course  of 
discharge.  The  comparison  between  limiting  cases  of  the  EDL  is 
made  using  equivalent  rate  constants  for  each  case.  It  remains  an 
open  issue  whether  different  rate  constants  could  be  selected  for 
the  Gouy-Chapman  model  that  would  provide  better  agreement 
with  experimental  data.  Whether  these  results  are  an  indication 
that  the  double  layer  structure  within  the  battery  studied  for  this 
work  tends  to  have  more  immobile  charge  in  the  Stern  layer  than 


from  theory  and  experimental  data.  Increasing  overpotential  is  modeled  as  a  linear 
function  of  depth  of  discharge. 


mobile  charge  in  the  diffuse  layer  is  left  as  an  interpretation  for  the 
reader.  The  Gouy-Chapman  results  also  diverge  from  the  experi¬ 
mental  data  for  the  4.8C  current  for  the  same  reasons  discussed  in 
the  context  of  the  Helmholtz  case. 


5.2.  Empirical  resistive  effects 

It  is  clear  from  the  preceding  results  that  some  effects  which 
may  be  important  for  accurately  predicting  cell  terminal  voltage 
during  high  discharge  currents  have  not  been  included  in  the 
theoretical  development  of  the  model.  As  discussed  previously,  the 
resistive  reactant  nature  of  iron  phosphate  results  in  increasing 
overpotential  with  respect  to  depth  of  discharge.  The  difference 
between  model  output  and  experimental  data  is  plotted  in  Fig.  10, 
for  the  range  of  0.8— 1.8  Ah  during  discharge  (roughly  the  middle 
third  of  the  discharge  curve). 

A  fit  of  the  overpotential  not  accounted  for  in  the  model  is  also 
shown  in  Fig.  10.  The  region  for  parameter  identification  is 
restricted  to  the  plotted  region  of  0.8  and  1.8  Ah.  This  is  because  the 
middle  portion  of  the  discharge  curve  is  less  prone  to  errors  from 
the  effects  of  incorrect  initial  conditions  that  may  cause  large  errors 
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in  predicted  voltage  at  the  beginning  and  end  of  the  discharge 
process.  These  potentially  large  errors  arise  because  the  open- 
circuit  voltage  varies  rapidly  with  respect  to  state  of  charge  in 
these  regions.  As  one  may  observe  from  Fig.  10,  the  unmodeled 
overpotential  is  linear  with  respect  to  charge  removed.  The  slope  of 
the  unmodeled  overpotential  is  of  use  as  a  fitting  parameter,  so  a 
linear  least  squares  problem  is  formulated  as 

Vrr  =  VrrQ  +  l/rrl  Q  (33) 

where  Vn  is  the  unmodeled  overpotential,  Vn$  is  the  initial 
unmodeled  overpotential  (the  value  at  Q  —  0.8  Ah),  and  VVi  is  the 
slope  of  the  overpotential  versus  charge.  When  solving  the  least 
squares  problem,  Vn  is  a  vector  of  values  corresponding  to  the 
difference  between  the  experimental  data  and  model  predicted 
voltage,  Vnfi  is  a  constant  that  is  subtracted  from  the  vector  Vn  to 
form  a  vector  Vrr,  and  Q.  is  a  vector  of  charge  values  corresponding 
to  Vn-.  The  solution  to  the  least  squares  problem  is  given  by 

Vn-,i;  =  VrrQ7  (oojy1  (34) 

where  the  slope  Vn-,1  is  in  terms  of  V  Ah-1.  Dividing  by  the 
magnitude  of  the  discharge  current  gives  a  slope  in  terms  of 
mO  Ah-1,  as  documented  in  Table  2. 

When  comparing  the  unmodeled  overpotential  on  a  voltage 
basis,  there  is  a  wide  discrepancy  between  the  various  discharge 
rates  reported  in  Table  2.  However  when  the  rate  of  voltage  change 
is  scaled  based  upon  the  magnitude  of  the  discharge  current,  the 
values  for  each  discharge  rate  show  much  less  variation.  Based  on 
this  analysis,  a  resistance  that  varies  with  respect  to  state  of  charge 
is  included  in  the  model  as 

Rn  =  6.3  Q  (35) 

where  the  value  of  the  overpotential  slope  from  the  4.8C  discharge 
has  been  used,  and  Qis  measured  from  an  initial  value  of  zero  at  the 
beginning  of  discharge  to  the  final  capacity  value  at  the  end  of 
discharge.  This  is  because  it  tends  to  have  the  largest  overpotential 
and  thus  the  greatest  signal  to  noise  ratio,  which  is  a  benefit  during 
the  parameter  identification  process.  The  overpotential  slopes  for 
the  C/3  and  C/1.2  are  more  likely  to  be  corrupted  by  the  presence  of 
noise  in  the  data. 

The  comparison  between  model  and  experiment  when  the 
resistive  reactant  effect  of  iron  phosphate  is  included  is  shown  in 
Fig.  11.  In  general  the  agreement  is  improved  more  for  the  higher 
rates  of  discharge,  since  the  resistance  has  greater  impact  on  the 
cell  voltage  at  higher  rates.  The  capacity  at  the  4.8C  rate  is  still 
underpredicted  due  to  neglect  of  thermal  effects.  Most  likely  in¬ 
ternal  heating  has  raised  the  cell  core  temperature  enough  to  cause 
an  elevated  diffusion  coefficient  and  enable  improved  utilization  of 
active  material. 


Fig.  11.  Comparison  of  model  output  (lines)  and  experimental  data  (symbols)  for  a 
range  of  current  rates  using  the  Helmholtz  model  with  resistive  reactant  effect 
included.  The  charge  values  along  the  x  axis  are  obtained  as  the  integral  of  the  constant 
current  rate  magnitude  with  respect  to  time.  Parameter  values  are  taken  from  Table  1. 


potential,  kinetic  overpotential,  and  solid  lithium  concentration  in 
each  electrode  has  been  presented.  The  primary  conclusion  is  that 
the  presence  of  diffuse  charge  in  the  EDL  causes  greater  initial 
overpotential  and  greater  variation  of  overpotential  with  respect  to 
battery  state  of  charge  when  modeled  in  an  intercalation  battery. 
We  make  this  conclusion  based  on  the  inclusion  of  aspects  of 
interfacial  physics  that  have  been  neglected  in  prior  work  and  the 
comparison  between  two  cases  of  the  EDL.  This  analysis  also  pro¬ 
vides  insight  for  understanding  the  Butler— Volmer  representation 
of  interfacial  charge  transfer  typically  applied  within  battery 
models.  Rather  than  interpreting  the  Butler-Volmer  law  as  a 
relationship  between  bulk  properties  of  the  solid  and  liquid,  it  may 
be  conceptualized  as  governing  the  relationship  between  the  re¬ 
action  rate  and  the  Stern  layer  voltage  if  diffuse  charge  is  not 
present. 

Furthermore,  modifications  have  been  proposed  to  the  theo¬ 
retical  current-voltage  relationship  of  the  EDL  and  electroneutral 
liquid.  It  has  been  shown  that  a  single  resistance  value  that  in¬ 
creases  with  depth  of  discharge  improves  the  agreement  with 
experimental  data  for  a  range  of  galvanostatic  discharge  experi¬ 
ments  from  C/3  to  4.8C.  The  model  presented  in  this  paper  captures 
essential  features  describing  capacity  such  as  the  amount  of 
cyclable  lithium,  active  material  volume  fraction,  and  rate  limita¬ 
tions.  Its  main  advantage  is  the  computational  simplicity  it  offers 
while  still  providing  agreement  with  galvanostatic  discharge  data. 
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Appendix 

Equilibrium  voltages  with  respect  to  a  Li/Li+  reference  must  be 
obtained  experimentally  for  the  anode  and  cathode  as  part  of  the 
usual  procedure  for  constructing  an  electrochemical  model.  The 
empirical  functions  from  Ref.  [25]  give  good  agreement  for  the  cells 
of  this  work.  For  the  anode, 
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U3„(a)  =  0.6379  + 0.5416  exp(-305.53  a) 
+  0.044  tan  h(-?-^s9g58) 


-  0.6875  tan  h 

-  0.0175  tan  h 


a+  0.0117\ 
0.0529  ) 
fa  -  0.5692\ 
\  0.0875  ) 


where  a  is  the  degree  of  lithiation  of  the  anode,  defined  as  the 
surface  concentration  divided  by  the  saturation  concentration.  For 
the  cathode, 


Ua(b)  =  3.4323  -  0.8428  exp  (  -  80.2493(1  -  b)1'3198) 

-  3.247  x  10  6  exp  (20.2645(1  -  b)3'8003)  +  3.2482 
x  10  6  exp  (20.2646(1  -  b)3'7995) 

(37) 


where  b  is  the  degree  of  lithiation  of  the  cathode,  defined  in  a 
fashion  similar  to  the  anode. 
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